Isobar of an ideal Bose gas within the grand canonical ensemble 
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We investigate the isobar of an ideal Bose gas confined in a cubic box within the grand canonical 
ensemble, for a large yet finite number of particles, N. After solving the equation of the spinodal 
curve, we derive precise formulae for the supercooling and the superheating temperatures which re- 
veal an iV -1 / 3 or A r_1//4 power correction to the known Bose-Einstein condensation temperature in 
the thermodynamic limit. Numerical computations confirm the accuracy of our analytical approxi- 
mation, and further show that the isobar zigzags on the temperature- volume plane if N > 14393. In 
particular, for the Avogadro's number of particles, the volume expands discretely about 10 5 times. 
Our results quantitatively agree with a previous study on the canonical ensemble within 0.1% error. 
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I. INTRODUCTION 

A classic paper by Anderson in 1972 goes with the ti- 
tle, More Is Different [l| , which characterizes the notion 
of 'emergence': the way complex systems and patterns 
arise out of a multiplicity of relatively simple interac- 
tions. One relevant question is then, How many is differ- 
ent? To answer the question, we may consult quantum 
statistical physics where the key quantity is the partition 
function. Once we know the exact expression of the par- 
tition function, we can compute various physical quan- 
tities. For example, when the partition function in the 
grand canonical ensemble, Z(T, V, z), depends on three 
variables (temperature, volume, fugacity), the pressure 
and the average number of the particles are given by 

P = k B Td v \nZ(T,V,z) , N = zd z In Z(T, V, z) , (1) 

where k B denotes the Boltzmann constant. If the system 
is finite, due to the analytic property of the partition 
function, the physical quantities which are given as a 
fraction between the partition function and its deriva- 
tives, like ([T]), cannot feature any mathematical singular- 
ities. On the contrary, infinite systems may do so. In 
this way, it seems that, More is the same; infinitely more 
is different @. 

Viewing the Avogadro's number, Na — 6 x 10 23 , as 
an enormous quantity might well suggest to take the in- 
finity limit or the thermodynamic limit: the limit of the 
large volume and the large number of particles with the 
density held fixed Q. Essentially due to the quantum 
commutation relation, [x,p] = ih, the reduced Planck's 
constant, h, is positioned inside the expression of the 
partition function along with the volume, V, gencrically 
through the combination, V/h 3 , where the power of h 
corresponds to the dimension of the space. This implies 
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that the large volume limit may be traded with the clas- 
sical limit H 0, and hence special care should be taken 
while considering the thermodynamic limit, in order to 
preserve any quantum nature [^-u|. Further, since taking 
the thermodynamic limit and taking the derivatives do 
not commute in general, desirably it is safer to take the 
thermodynamic limit only at the end of computation. 

Recently two of the authors investigated the isobar of 
an ideal Bose gas confined in a box within the canonical 
ensemble, without assuming the thermodynamic limit 0] . 
Numerical computations based on the exact expression of 
the corresponding canonical partition function revealed 
that, if the number of particles is equal to or greater 
than a certain critical value, which turns out to be 7616 
for the 'cubic' box, the isobar zigzags featuring '5-shapc' 
on the (T, F)-plane (c/. FIG. [2] in the present paper). 
The two turning points on the 5-shaped isobar are nat- 
urally identified as the 'supercooling' (T*,V*) and the 
'superheating' (x**^**) points. Between the supercool- 
ing and the superheating temperatures, T* < T < T**, 
the volume becomes triple-valued. Since all the physi- 
cal quantities are functions of the temperature and the 
volume, every physical quantity itself is triple- valued be- 
tween the two temperatures and changes discontinuously 
on isobars as the temperature increases. In fact, any tem- 
perature derivative restricted on isobars diverges at the 
points with the universal singularity exponent, 1/2 [9[. 
In this way imposing the 'constant pressure constraint,' 
a discrete phase transition was for the first time realized 
in a finite system, derived ab initio from the correspond- 
ing partition function. 

However, due to the limitation in our computational 
power (supercomputer, SUN B6048), the numerical anal- 
yses performed in Refs.@, Q were restricted to the parti- 
cle numbers not greater than one million. In particular, 
the separation between the supercooling and the super- 
heating temperatures gets wider as the number of parti- 
cles increases within the range, 7616 < N < 10 6 . Hence, 
it was not clear, what would happen for much larger num- 
ber of particles, or closer to the thermodynamic limit. 
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It is the purpose of the present paper, first to verify 
the same feature of the ideal Bose gas within the grand 
canonical ensemble, both analytically and numerically; 
and second to address rigorously its thermodynamic 
limiting behavior. 

Basically we set to analyze the following equation [l(3| 
which shall be derived from the grand canonical partition 
function of the ideal Bose gas: 



dT 

av 



= 0. 



(2) 



P,N 



This condition is equivalent to the usual definition of the 
spinodal curve fl(H13l |. 



dP 

av 



= 0. 



(3) 



T,N 



and must be met at the supercooling and the superheat- 
ing points on isobars. 



II. ANALYSIS 

Essentially due to the non-relativistic dispersion rela- 
tion, E — p 2 1 (2m), where m is the mass of the particle, 
the grand canonical partition function of the ideal Bose 
gas confined in a cubic box is essentially a two-variable 
function depending on the fugacity, z, and the combina- 
tion of temperature and volume, TV 2 / 3 . Specifically we 
set, as for the two fundamental variables in our analysis, 



2m kv 



(TV 



2/3) 



: = — In ; 



(4) 



In terms of these, the grand canonical partition function 
reads 



\nZ{e,<r) = hill- 



-en —a 



(5) 



With the Dirichlet boundary condition which we deliber- 
ately impose, n = (m, 712,713) <E N" 3 is a positive integer- 
valued lattice vector, such that the lowest value of n 2 is 
the spatial dimension, 3, and a is bounded from below 



cr > — 3e , 



(0) 



while e is positive. Searching for spinodal curves near to 
the thermodynamic limit, we shall be interested in the 
small e region. 

It is useful to note, for the computation of various phys- 
ical quantities such as ([I]), 

Td T \ v , z = lVdv\ TtZ = -Ed e , zd z \ TV = -d„ . (7) 

It follows that the number of particles ([T]) reads 

N(e, a) = -cUn 2(6,(7), (8) 



and the formula of the pressure ([I]) is equivalent to 
T P (e,a) :=(^)h s TP-i = In Z( £ , a 

(9) 

Being a combination of T and P, this dimensionless quan- 
tity, T P , can determine the physical temperature on an 
arbitrarily given isobar. Similarly we may define a di- 
mensionless "volume" , 



V P (e,a) := {^PY V = [-ld s lnZ(e,a)] 
and another dimensionless "temperature" , 



T p {e,a) 



2m 7 rp (V\S 



-e ^dcr lnZ(e, a) 



(10) 



(11) 



As we already wrote, N, T P , V P and T p arc functions of 
the two variables, e, a only. They satisfy identities, 



eT P (e ) a) = [V P (s,a)Y 



eT P (e,a) = [N(e,a)]-* 



(12) 

Now the spinodal curve ^ is positioned on the (e, em- 
plane to satisfy diV(e,cr) = and dT P (e,cr)=0, such 
that the following linear equation must admit a nontrivial 
solution, 

/ 





-( 







d £ d a lnZ 



if- In Z 



[le^d E + d 2 )\nZ d E d a hiZ 



de 
da 



(13) 

It follows that the 2x2 matrix in (flU)) must be singular, 

d e d a lnZ <9 2 lnZ 
^ds + d^lnZ d E d a \nZ 



$ := det 



0. (14) 



This algebraic equation determines the spinodal curve on 
the [e, cr)-planc. Further, it is straightforward to show 

that the determinant is proportional to QO 



N 



d In Tp 
d In V P 



N 



3(9J lnZ) 2 Var(n 2 



X $ . 



where Var(7i 2 ) is our shorthand notation for 
Var(n ) ■ 



/ d B d a \nZ \' 



(15) 



(16) 



which can be identified as the variance of ft 2 with 
respect to the probability distribution proportional to 
sinh _2 (ien 2 + \a) [ll|. Hence, Var(7i 2 ) is positive 
definite and the vanishing of the determinant is, as 



expected, equivalent to the vanishing of 



Oi 



N 



main task is to solve (|14[) and express the solutions in 
terms of the more physical variables, N, T P , V P , Tp 
using ([8|), ([9|), ([TO]) . (fTTjl . Our numerical solutions are 
depicted in FIGQ] and FIG|H along with an analytic 
approximation which we discuss below. 
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FIG. 1: The supercooling and the superheating spinodal 
curves on the (TV -1 , 7p/7]? BC )-plane (lower and upper 
curves respectively). The dotted curves are from the numer- 
ical computations based on the exact formulae ([8]), (JHJ), fl4)l . 
The solid lines correspond to our analytic approximation (|34[1 , 
(|35f) for large TV. A pair of spinodal curves start to develop at 
TV = TV C ~ 14392.4 (iV c _1/4 ~ 0.0912991) which is comparable 
to the critical number of the canonical ensemble, 7616 Q. 
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FIG. 2: Isobar curves on the (Tp/T p bec , V P /TV)-plane. They 
zigzag featuring '5-shape' if 14393 < TV < oo. 



Analytic approximation 

Our analytic analysis starts with the following expres- 
sion for the derivatives of the partition function, 



(d e y(d„) k - l lnZ(e,a) 



t2\i (-irc fc „ 



(17) 



Here k > max(Z,l), and Cfc j(1 are positive integers that 
are determined by a recurrence relation, 



Ck+l.a — aCk,a + (a — l)Cfc, a _i , 



(18) 



with the initial value, Ci.i = 1. The recurrence relation 
comes from the expansion, 

HI -e-n^-ELiT^try- (19) 

Taking an ^-derivative of the right hand side of the equal- 
ity in (fT^|) leads to (fT5|) . Further, it is useful to note 



d 2 lnZ( £ ,a) = E 



riGF 



(e»1 2 +<x) 2 



£(i)*coslT 



2 / en 2 +er 



fc=2 



This expression is due to an identity, 



(20) 



sinh 2 (a;) 

= (i)"sinh- 2 (f ) -£?=i(^co S h- 2 (§ ) (21) 
= -- 2 -E,°:i(4) J cosh- 2 (f), 

which holds for an arbitrary positive integer, p. Taking 
p to infinity gives the second equality in (1211) . 



In order to compute the sums in (j 1 T[> . we adopt the 
following scheme of analytic approximation: 



1. Introduce a cut-off, A > 3, for the lattice sum, 

E ^ 2 ) = E /( £ " 2 )+ E /( e " 2 ) • 



(22) 



n 2 <A 



n 2 >A 



2. Approximate the last term by an integral, 

/•OO 

E /M 2 ) ^ E /( £ " 2 )+ / d:E (f^ 4 ^ - TrO/M- 



rleN 3 



n 2 <A 



(23) 



3. Piti (7 = — e/x wif/i a new variable, fx. From 
fi < 3. 

4. Keep only the dominant singular terms in the power 
series expansion of ()23|) in e, which are manifestly 
cut-off independent. Allow /z to be expandable in e 
with an arbitrary leading power. 



The approximation (|23|) can be traced back to an identity, 



+IE„ ez /(^ 2 )]-|/(o), 

(24) 

where the first two sums on the right hand side of the 
equality can be approximated by integrals in three or two 
dimensional spherical coordinates, and the remaining 
part may be neglected for small e (see 0| and references 
therein). 
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With the constants, 

a s := / °° dx^ = T(s + l)C(s + 1) , 
b := / °° dx y/x cosh -2 (a;) , 

and the estimations 3 , 

Q dx W=T ~ JeA dx I^=TF - - ln£ ■ 

our scheme enables us to compute 



(25) 



(26) 



cases, h = 1 and h = 2, admit solutions, and hence there 
are two spinodal curves as follows. 



• On the (e, ^)-plane. 
* constant /i ~ /i* line with ft = 1, satisfying 



^nSN 3 (n 2 -^*) 2 

Numerically we get 

\i* ~ 2.61873. 



(31) 
(32) 
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>(3-M) 2/ -t 



03£ 2 

2 

3-7T 



£ In £ , 



OiE 2 + 4f £ ln£ , 



(3^F 



>3 (n 2 -^) 2 



r -2 _ (2+V2)7r & , 



(27) 

where (g(e))- n denotes a part of the series expansion 
of g(e) in e which is at least (— n)-th order singular, for 
example, 



_ 3 

£ 2 



£ -2 + £~ X + 1 + S)_, 

2 

£ _1 + l + e)-i = e~i 



I -2. 

(£ 2 



(28) 



etc. 



Especially for <9 2 lnZ, it is important to note that the 
sum, X)n 2 >3("- 2 — > converges, since 



n 2 >A (n 2 — u) 2 — Z-^n 



>A (n 2 -|H) 2 

3tt 



fOO , -fa; 2 

Ja dx 



The numerical values of the constants are 



(29) 



oi = ^C(f) - 2.31516 , ai = ~ 1.64493 , 
a| = ^C(f) ^ 1-78329, 6 - 0.758128. 

(30) 

Having the expressions (f27|) , we now proceed to solve the 
spinodal curve condition (|14[). Since the indices, n of the 
symbol, ( • }_„ appearing in (|27p are various, letting the 
leading singular term of be order of £~ ft , wc need to 
separately consider the following nine possible cases: 



** Linear line with h = 2, 

> In terms of the physical variables, TV, T P , V P , T p 
* Supercooling spinodal curve, for ft = 1, 



(33) 



1 + -2- 

^60 



(7 -BEC)5 /T BEC 



r;/7 -BEC 

v* ~ (r p BEC /r P BEC ) 1 (jv + f r p BEC jvf in 
r;/r p BEC ~ 1 + f t p bec in jv. 



AT-3 

AT 



(34) 



** Superheating spinodal curve, for ft = 2, 



vr ~ 8 



1 



15 \ 4 /fBEC 



(r P BEC ) -5 ivt , 



(35) 



In the above, T P EC and "7p BEC denote two constants, 

T P BEC = [C(f)] _ * ^ 1.02781 , 

V EC = i [C(§)]~ f - 0.671253, 



(36) 



which correspond to the well-known Bosc-Einstein con- 
densation temperatures for the variables, T P j9j and T P 
(fTTj) , the definitions of which we recall here, 



Tp := 



-2lSl.\ 5 b TP~ 

^2fi2 1 tB- 1 ^ 



J r, 2 
VV 



-T- 2m u rp I V_ \ 3 



2m 

TT 2 K 2 



P) 5 ^ = N 



(37) 



III. DISCUSSION 



ft < 1 , ft = 1 , 1< ft < | 



ft = & 

" 4 ' 4 



5 < /, < 3 
, v. n „ , 



<ft<2, ft = 2, 2<ft. 



Keeping only the two dominant terms in (|27|) for each 
case, it is straightforward to check that only the two 



As computable from our analytic expressions, (|34p and 
(|35|) . the separation between the supercooling and the su- 
perheating temperatures becomes maximal, if the num- 
ber of particle is equal to 



N MAX - j^pp [((§)] 



2.32890 x 10 b 



(38) 
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This also agrees with the numerical result in FIG. [T] as 
(TVmax)~ 1/4 ~ 0.0255984. When the number of particles 
exceeds this critical value, the two temperatures, T P and 
T**, - satisfying 7p BEC < T* < T** - gets closer, and 
eventually converges to the BEC temperature, Tp EC ([36f , 
in the thermodynamic limit. That is to say, TV MAX is the 
critical number for the thermodynamic limit to work. 
The ratio of the two volumes, 

V;/V* p * ~ (^)* [a^V'Ni ~ 0.118511 x Ni , (39) 

enables us to estimate the discrete volume expansion rate 
at the liquid-gas type phase transition. For the Avo- 
gadro's number, Na — 6.022 1 4 x 10 23 , the volume ex- 
pansion rate fl35J) gives V*/V** ~ 104399. Thus, the 
ideal Bose gas made up of the Avogadro's number of 
particles expands its volume discretely about 10 5 times 
during the phase transition. This is a genuine finite effect 
of the Avogadro's number, which cannot be seen directly 
in the thermodynamic limit where V P /V%* — > oo. 

Our numerical computations based on the exact for- 
mulae quantitatively agree with the canonical ensemble 
results [|] for TV = 10 5 and TV = 10 6 within 0.1% error, 
though the minimum (natural) numbers required for the 
emergence of the spinodal curves are different, 14393 vs. 
7616. 





Grand canonical 


Canonical 


TV = 10 5 


(1.041,1.043) 


(1.0410, 1.0424) 


TV = 10 6 


(1.0348, 1.0364) 


(1.034, 1.036) 



TABLE I: Quantitative agreement between the canonical and 
the grand canonical results, within 0.1% error. 



In this work, we have focused on the Dirichlct 
boundary condition. Alternatively imposing periodic 
or Neumann boundary condition brings out a volume 
independent ground state energy which, as shown in Q, 



causes a thermodynamic instability at low temperature 
near absolute zero (see also [l4|). This further implies 
that, under the alternative boundary conditions, periodic 
or Neumann, the isobar on the (T Pl V P /N) plane is of 
'C-shape', rather than of the zigzagging '5-shape' as in 
FIGlU Namely there is a nontrivial lower bound in T P 
of the isobar, above which the volume is always double- 
valued. In the thermodynamic limit, the lower bound 
converges to T^ EC , and the isobar eventually becomes 
independent of the boundary conditions, identical to the 
case of TV = oo in FIGH except for V P /N = 0. When 
Vp/TV = 0, under the periodic or Neumann boundary 
condition, T P may assume any value which is greater 
than or equal to T P EC (as anticipated in Fig. 12.8 of 0]), 
while under the Dirichlet boundary condition, it is quite 
the opposite, < T P < JJ EC \ as depicted in FIGH 

In conclusion, we have shown, both numerically 
and analytically, that the isobar of the ideal Bose gas 
zigzags on the temperature-volume plane, qualitatively 
featuring the liquid-gas transition, if TV > 14393. This is 
an emergent phenomenon of the finitely many bosonic 
identical particles. We have derived the precise formulae 
for the the two turning points: supercooling (|34|) and 
superheating (|35[) . Our formulae reveal an TV -1 / 3 or 
TV -1 / 4 power correction to the BEC temperature and 
enable us to estimate the volume expansion rate, (1551) . 
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